voronoi/gen_voronoi_3D.py (2)

20-22: Avoid mutating the caller’s input array.

points = initial_points aliases the input, so the relaxation mutates initial_points and the “initial points” plot becomes incorrect. Copy the array before iterating.

🛠️ Suggested fix

-    points = initial_points
+    points = np.array(initial_points, copy=True)

---

61-63: Create the output directory recursively.

os.mkdir fails if ../runs doesn’t exist. Use os.makedirs(..., exist_ok=True) instead.

🛠️ Suggested fix

-    if os.path.exists(output_dir) == False:
-        os.mkdir(output_dir)
+    os.makedirs(output_dir, exist_ok=True)

examples/3D_random_sphere_array/case.py (2)

1-3: Unused import: math

The math module is imported but never used. All mathematical operations use numpy instead.

Suggested fix

 import json
-import math
 import numpy as np

---

36-38: Consider adding basic validation for the sphere location file.

If sphere_array_locations.txt is missing or has fewer than 3 columns, the script will fail with an unclear error when accessing sphere_loc[i, 0:3] in the loop below. A brief check could improve usability:

Suggested validation

 # load initial sphere locations
-sphere_loc = np.loadtxt("sphere_array_locations.txt")
+sphere_loc = np.loadtxt("sphere_array_locations.txt", ndmin=2)
+if sphere_loc.ndim != 2 or sphere_loc.shape[1] < 3:
+    raise ValueError("sphere_array_locations.txt must contain at least 3 columns (x, y, z)")
 N_sphere = len(sphere_loc)

src/simulation/m_data_output.fpp (3)

1018-1054: Avoid writing the beta slot twice in the filtered MPI loop.
When q_filtered_wrt is true and beta is present, the loop writes sys_size+1 and then the “Correction” block writes it again, risking overlap/corruption.

🐛 Suggested fix

            else if (q_filtered_wrt .and. (t_step == 0 .or. t_step == t_step_stop)) then
                do i = 1, alt_sys
+                    if (present(beta) .and. i == sys_size + 1) cycle
                     var_MOK = int(i, MPI_OFFSET_KIND)
                     disp = m_MOK*max(MOK, n_MOK)*max(MOK, p_MOK)*WP_MOK*(var_MOK - 1)
                     call MPI_FILE_SET_VIEW(ifile, disp, mpi_p, MPI_IO_DATA%view(i), &
                                            'native', mpi_info_int, ierr)
                     call MPI_FILE_WRITE_ALL(ifile, MPI_IO_DATA%var(i)%sf, data_size*mpi_io_type, &
                                             mpi_io_p, status, ierr)
                end do

---

117-123: Serial path drops filtered outputs.
When parallel_io is false, the filtered/statistics arguments are ignored, so q_filtered_wrt cannot produce outputs in serial runs. Please extend the serial path or explicitly reject this configuration.

💡 Suggested fix (also update s_write_serial_data_files signature)

-        if (.not. parallel_io) then
-            call s_write_serial_data_files(q_cons_vf, q_T_sf, q_prim_vf, t_step, bc_type, beta)
-        else
+        if (.not. parallel_io) then
+            if (q_filtered_wrt .and. present(filtered_fluid_indicator_function)) then
+                call s_write_serial_data_files(q_cons_vf, q_T_sf, q_prim_vf, t_step, bc_type, beta, &
+                                               filtered_fluid_indicator_function, &
+                                               stat_q_cons_filtered, stat_filtered_pressure, &
+                                               stat_reynolds_stress, stat_eff_visc, stat_int_mom_exch)
+            else
+                call s_write_serial_data_files(q_cons_vf, q_T_sf, q_prim_vf, t_step, bc_type, beta)
+            end if
+        else
             call s_write_parallel_data_files(q_cons_vf, t_step, bc_type, beta, &
                                              filtered_fluid_indicator_function, &
                                              stat_q_cons_filtered, stat_filtered_pressure, &
                                              stat_reynolds_stress, stat_eff_visc, stat_int_mom_exch)
         end if

---

957-963: Filtered MPI buffers can remain uninitialized.
s_initialize_mpi_data_filtered is only called under ib, so q_filtered_wrt with ib = .false. (or file-per-process) leaves filtered views/buffers uninitialized. Move the call outside the ib branch or assert unsupported combinations.

src/simulation/m_volume_filtering.fpp (3)

585-594: Finite-difference stencil can read out of bounds.
The loop runs i=0..m, j=0..n, k=0..p while indexing i+l with l=-4..4, which can access beyond bounds near edges. Restrict loop bounds or use buffered indices.

🐛 Suggested fix (bounds-safe loop)

-        $:GPU_PARALLEL_LOOP(collapse=3, copyin='[fd_coeffs]')
-        do i = 0, m
-            do j = 0, n
-                do k = 0, p
+        $:GPU_PARALLEL_LOOP(collapse=3, copyin='[fd_coeffs]')
+        do i = 4, m-4
+            do j = 4, n-4
+                do k = 4, p-4
                     $:GPU_LOOP(parallelism='[seq]')
                     do l = -4, 4
                         grad_fluid_indicator(1)%sf(i, j, k) = grad_fluid_indicator(1)%sf(i, j, k) + fd_coeffs(l)*fluid_indicator_function%sf(i + l, j, k)/dx(i)
                         grad_fluid_indicator(2)%sf(i, j, k) = grad_fluid_indicator(2)%sf(i, j, k) + fd_coeffs(l)*fluid_indicator_function%sf(i, j + l, k)/dy(j)
                         grad_fluid_indicator(3)%sf(i, j, k) = grad_fluid_indicator(3)%sf(i, j, k) + fd_coeffs(l)*fluid_indicator_function%sf(i, j, k + l)/dz(k)
                     end do
                 end do
             end do
         end do

---

785-786: Guard against zero filtered density.
q_cons_filtered(1)%sf can be zero inside particles, causing division-by-zero/NaNs in Reynolds stress.

🐛 Suggested fix

+        real(wp), parameter :: eps_rho = 1.0e-12_wp
...
-                            reynolds_stress(l, q)%sf(i, j, k) = reynolds_stress(l, q)%sf(i, j, k) &
-                                                                - (q_cons_filtered(momxb - 1 + l)%sf(i, j, k)*q_cons_filtered(momxb - 1 + q)%sf(i, j, k)/q_cons_filtered(1)%sf(i, j, k))
+                            reynolds_stress(l, q)%sf(i, j, k) = reynolds_stress(l, q)%sf(i, j, k) &
+                                                                - (q_cons_filtered(momxb - 1 + l)%sf(i, j, k)*q_cons_filtered(momxb - 1 + q)%sf(i, j, k) / &
+                                                                   max(q_cons_filtered(1)%sf(i, j, k), eps_rho))

---

1553-1560: Filter back-transform divides by near-zero fluid indicator.
Several assignments divide by filtered_fluid_indicator_function%sf, which can be tiny near/inside particles. Add an epsilon floor (or conditional zeroing) for stability.

🐛 Suggested fix (apply to all similar divisions)

+        real(wp), parameter :: eps_phi = 1.0e-12_wp
+        real(wp) :: denom
...
-                        q_cons_filtered(l)%sf(i, j, k) = data_real_3D_slabz(i + 1, j + 1, k + 1)*fft_norm/filtered_fluid_indicator_function%sf(i, j, k)
+                        denom = max(filtered_fluid_indicator_function%sf(i, j, k), eps_phi)
+                        q_cons_filtered(l)%sf(i, j, k) = data_real_3D_slabz(i + 1, j + 1, k + 1)*fft_norm/denom

Also applies to: 1602-1607, 1652-1656, 1704-1708

src/simulation/p_main.fpp (1)

59-64: Ensure volume-filtering resources are finalized.
Initialization is added here, but the matching finalization calls are not visible. Please confirm these are invoked in s_finalize_modules() when volume_filter_momentum_eqn is enabled.

#!/bin/bash
# Verify volume-filtering finalizers are called in shutdown paths.
rg -n "s_finalize_.*filter|s_finalize_fftw_explicit_filter_module" src/simulation
rg -n "volume_filter_momentum_eqn" src/simulation/m_start_up.fpp src/simulation/p_main.fpp

src/common/m_ib_patches.fpp (1)

617-633: Avoid shared module state inside the GPU loop (thread-safety).

s_convert_cylindrical_to_cartesian_coord writes module-level cart_y/cart_z inside a parallel region. This is the same race previously flagged; please make the conversion thread‑local (e.g., using f_convert_cyl_to_cart and local temps) instead of module variables.

src/common/m_mpi_common.fpp (1)

1604-1613: Slab decomposition validation relies on generic error message.

The validation for slab decomposition divisibility is present via the ierr flag check at lines 1731-1734, but the error message is generic ("Unsupported combination of values..."). Consider adding a more specific error message for the slab decomposition case to help users diagnose the issue. This is a duplicate of a past review suggestion.

src/simulation/m_data_output.fpp (2)

85-115: Reduce argument count in s_write_data_files.
This signature now exceeds the ≤6-argument guideline; consider a derived-type params struct for the filtered/statistics payload. As per coding guidelines, ...

---

801-814: Reduce argument count in s_write_parallel_data_files.
This signature exceeds the ≤6-argument guideline; bundle the optional filtered/statistics fields into a derived-type params struct. As per coding guidelines, ...

src/simulation/m_volume_filtering.fpp (1)

721-729: Trim parameter lists for filtering routines.
s_setup_terms_filtering and s_filter_batch exceed the ≤6-argument guideline; wrap these inputs into a derived-type params object for maintainability. As per coding guidelines, ...

Also applies to: 1037-1046

src/post_process/m_data_input.f90 (1)

664-669: Reduce argument count via a params struct.

s_read_parallel_filtered_data takes 8 arguments, exceeding the project guideline (≤6). Consider a small derived type (e.g., type(mpi_io_params)) to pass the MPI offsets as a bundle.

As per coding guidelines, keep routine arguments ≤ 6 by grouping related parameters.

toolchain/mfc/run/case_dicts.py (1)

118-143: Large parameter dictionary expansion for IB patches.

Expanding ib_id range from 10 to 1000 generates thousands of dictionary entries. While this aligns with the PR objective to support up to 1000 IB patches, consider the impact on:

1. Schema validation performance (fastjsonschema compilation)
2. Memory footprint of the parameter dictionaries

This may be acceptable for the use case, but worth noting for future optimization if startup time becomes an issue.

src/simulation/m_global_parameters.fpp (1)

1048-1050: Consider documenting the volume_filter_size formula.

The formula sys_size + 1 + 4*(2*num_dims**2 + num_dims + E_idx + 1) accounts for multiple filtered quantities with 4 statistical moments each, but the calculation is not self-documenting. Consider adding a comment explaining the breakdown:

• sys_size: base conservative variables
• +1: filtered fluid indicator function
• 4 * (...): 4 statistical moments for each filtered quantity

src/simulation/p_main.fpp (1)

25-64: Add finalization for volume-filtering resources.

The new initializers are invoked under volume_filter_momentum_eqn, but there’s no corresponding cleanup in the finalize path. Please add finalizers (or confirm s_finalize_modules already handles them) to avoid leaking large allocations.

src/common/m_mpi_common.fpp (1)

341-389: Fix alt_sys so filtered MPI views cover the full variable range.

alt_sys is set to volume_filter_size, so the view-creation loop can miss appended variables. It should be sys_size + volume_filter_size and the loop should start at sys_size + 1.

🐛 Proposed fix

-        alt_sys = volume_filter_size
+        alt_sys = sys_size + volume_filter_size
...
-        do i = sys_size, alt_sys
+        do i = sys_size + 1, alt_sys

src/post_process/m_data_input.f90 (1)

679-709: Fix alt_sys so MPI filtered reads include base + filtered variables.

alt_sys is set to volume_filter_size, so NVARS_MOK and the read loop undercount when base variables are present. It should be sys_size + volume_filter_size.

🐛 Proposed fix

-        alt_sys = volume_filter_size ! filtered indicator, stats of: R_u, R_mu, F_imet, q_cons_filtered, pressure
+        alt_sys = sys_size + volume_filter_size ! filtered indicator + stats appended to base vars

src/simulation/m_ibm.fpp (3)

620-630: Fix off-by-one in image-point search bounds.
index+1 is accessed, so allowing index==m+buff_size can read past the coordinate arrays.

🐛 Proposed fix

-                    bound = m + buff_size
+                    bound = m + buff_size - 1
...
-                    bound = n + buff_size
+                    bound = n + buff_size - 1
...
-                    bound = p + buff_size
+                    bound = p + buff_size - 1

---

648-651: Use s_mpi_abort instead of error stop.
error stop can hang MPI ranks; switch to the MPI-aware abort. As per coding guidelines, use s_mpi_abort for error termination.

🐛 Suggested fix

-                            print *, "Increase buff_size further in m_helper_basic (currently set to a minimum of 10)"
-                            error stop "Increase buff_size"
+                            call s_mpi_abort("Increase buff_size further in m_helper_basic (currently set to a minimum of 10)")

---

1316-1320: Deallocate levelset fields even when store_levelset=F.
They are allocated unconditionally in initialization, so conditional deallocation leaks memory. As per coding guidelines, use the deallocation macro.

🐛 Proposed fix

         @:DEALLOCATE(ib_markers%sf)
-        if (store_levelset) then
-            @:DEALLOCATE(levelset%sf)
-            @:DEALLOCATE(levelset_norm%sf)
-        end if
+        @:DEALLOCATE(levelset%sf)
+        @:DEALLOCATE(levelset_norm%sf)

src/post_process/m_global_parameters.fpp (2)

123-124: Ensure volume_filter_size is initialized for all q_filtered_wrt cases.
It’s computed only under model_eqns == 2, but used for allocation whenever q_filtered_wrt is true; if q_filtered_wrt is enabled for another model, the size is undefined. Consider computing it after sys_size is finalized for all models, or validate that q_filtered_wrt is only used with model_eqns=2.

♻️ Suggested placement (after sys_size is finalized)

-            if (q_filtered_wrt) then
-                volume_filter_size = sys_size + 1 + 4*(2*num_dims**2 + num_dims + E_idx + 1)
-            end if
...
+        if (q_filtered_wrt) then
+            volume_filter_size = sys_size + 1 + 4*(2*num_dims**2 + num_dims + E_idx + 1)
+        end if

Also applies to: 698-700, 865-868

---

1058-1062: Deallocate extra MPI_IO_DATA entries for filtered outputs.
Nulling pointers alone leaks the allocated sf buffers. As per coding guidelines, deallocate with the macro when associated.

🐛 Proposed fix

             if (q_filtered_wrt) then
                 do i = sys_size + 1, volume_filter_size
-                    MPI_IO_DATA%var(i)%sf => null()
+                    if (associated(MPI_IO_DATA%var(i)%sf)) then
+                        @:DEALLOCATE(MPI_IO_DATA%var(i)%sf)
+                    end if
                 end do
             end if

src/simulation/m_data_output.fpp (3)

117-123: Serial path drops filtered outputs when parallel_io is false.

q_filtered_wrt data never reaches s_write_serial_data_files, so filtered outputs are silently skipped in serial runs. This matches a prior finding.

🔧 Suggested fix

-        if (.not. parallel_io) then
-            call s_write_serial_data_files(q_cons_vf, q_T_sf, q_prim_vf, t_step, bc_type, beta)
+        if (.not. parallel_io) then
+            if (q_filtered_wrt .and. present(filtered_fluid_indicator_function)) then
+                call s_write_serial_data_files(q_cons_vf, q_T_sf, q_prim_vf, t_step, bc_type, beta, &
+                                               filtered_fluid_indicator_function, stat_q_cons_filtered, stat_filtered_pressure, &
+                                               stat_reynolds_stress, stat_eff_visc, stat_int_mom_exch)
+            else
+                call s_write_serial_data_files(q_cons_vf, q_T_sf, q_prim_vf, t_step, bc_type, beta)
+            end if
         else
             call s_write_parallel_data_files(q_cons_vf, t_step, bc_type, beta, &
                                              filtered_fluid_indicator_function, &
                                              stat_q_cons_filtered, stat_filtered_pressure, &
                                              stat_reynolds_stress, stat_eff_visc, stat_int_mom_exch)
         end if

---

957-963: Filtered MPI buffers are still only initialized inside the ib branch.

When q_filtered_wrt is true but ib is false (or when file_per_process is true), the filtered MPI buffers/views remain uninitialized. This mirrors a prior finding.

---

1018-1054: Potential duplicate write of beta slot when q_filtered_wrt is enabled.

The q_filtered_wrt loop writes indices 1..alt_sys, which includes sys_size+1 when present(beta) sets alt_sys = sys_size + 1, and then beta is written again in the correction block. This is a previously reported overlap risk.

src/simulation/m_additional_forcing.fpp (2)

111-115: Momentum-energy term still assumes 3 components.

momxb, momxb+1, momxb+2 hard-code a 3D layout and break 1D/2D or nonstandard layouts. This matches a prior review note.

---

166-168: Energy forcing denominator is single-fluid and unguarded.

Division by q_cons_vf(contxb) can be zero and ignores multi-fluid mixture density, consistent with earlier review findings.

src/simulation/m_global_parameters.fpp (1)

1048-1050: volume_filter_size only set in model_eqns==2.

If q_filtered_wrt is enabled under other equation sets, volume_filter_size can be uninitialized and lead to invalid allocations. This repeats a prior review finding.

src/simulation/m_compute_statistics.fpp (1)

239-246: Moments divide by zero when ns <= 3 or variance is zero.

The formulas can divide by ns-1/2/3 and by Msn(2) when the variance is zero. This matches a prior review finding.

src/simulation/m_volume_filtering.fpp (3)

585-594: Stencil access can go out of bounds at domain edges.
i+l, j+l, k+l will underflow/overflow near boundaries unless halos are ≥4. Restrict loop ranges or guard indices.

---

785-786: Guard division by filtered density.
q_cons_filtered(1)%sf(i,j,k) can be ~0 in solid regions, causing NaNs. Add a small floor.

🐛 Suggested pattern

-                            reynolds_stress(l, q)%sf(i, j, k) = reynolds_stress(l, q)%sf(i, j, k) &
-                                                                - (q_cons_filtered(momxb - 1 + l)%sf(i, j, k)*q_cons_filtered(momxb - 1 + q)%sf(i, j, k)/q_cons_filtered(1)%sf(i, j, k))
+                            reynolds_stress(l, q)%sf(i, j, k) = reynolds_stress(l, q)%sf(i, j, k) &
+                                                                - (q_cons_filtered(momxb - 1 + l)%sf(i, j, k)*q_cons_filtered(momxb - 1 + q)%sf(i, j, k) &
+                                                                   / max(q_cons_filtered(1)%sf(i, j, k), eps_rho))

---

1553-1557: Protect divisions by filtered_fluid_indicator_function.
The denominator can be 0 in solid cells, leading to NaNs for q_cons_filtered, filtered_pressure, reynolds_stress, and eff_visc. Add an epsilon or conditional guard.

🐛 Suggested pattern (apply to all four sites)

-                        q_cons_filtered(l)%sf(i, j, k) = data_real_3D_slabz(i + 1, j + 1, k + 1)*fft_norm/filtered_fluid_indicator_function%sf(i, j, k)
+                        q_cons_filtered(l)%sf(i, j, k) = data_real_3D_slabz(i + 1, j + 1, k + 1)*fft_norm / &
+                            max(filtered_fluid_indicator_function%sf(i, j, k), eps_phi)

Also applies to: 1602-1606, 1652-1656, 1704-1708

toolchain/mfc/case_validator.py (1)

1104-1155: Reduce local variable count in check_periodic_forcing to satisfy pylint.
A small table-driven validation can drop several locals while keeping the same error messages.

♻️ Possible refactor sketch

-        u_inf_ref = self.get('u_inf_ref')
-        rho_inf_ref = self.get('rho_inf_ref')
-        P_inf_ref = self.get('P_inf_ref')
-        forcing_window = self.get('forcing_window')
-        forcing_dt = self.get('forcing_dt')
+        required = {
+            'u_inf_ref': "u_inf_ref (desired bulk velocity) must be specified for periodic_forcing",
+            'rho_inf_ref': "rho_inf_ref (desired bulk density) must be specified for periodic_forcing",
+            'P_inf_ref': "P_inf_ref (desired pressure based on bulk internal energy) must be specified for periodic_forcing",
+            'forcing_window': "forcing_window must be specified for periodic_forcing",
+            'forcing_dt': "forcing_dt must be specified for periodic_forcing",
+        }
+        for key, msg in required.items():
+            self.prohibit(self.get(key) is None, msg)

src/simulation/m_ibm.fpp (1)

154-228: Add explicit private lists to periodic-buffer GPU loops.
If the macro doesn’t implicitly privatize j/k/l, these loops can race. Please verify and add private='[j,k,l]' (apply to all six loops). As per coding guidelines, declare loop-local variables in GPU parallel macros.

♻️ Example (apply to all six loops)

-                $:GPU_PARALLEL_LOOP(collapse=3)
+                $:GPU_PARALLEL_LOOP(collapse=3, private='[j,k,l]')

src/simulation/m_volume_filtering.fpp (2)

721-729: Consider grouping arguments to keep subroutines ≤6 params.
Both s_setup_terms_filtering and s_filter_batch exceed 6 arguments. A derived-type params struct would keep signatures manageable. As per coding guidelines, consider consolidating related fields.

Also applies to: 1037-1045

---

92-95: Verify mixed‑precision slab buffers are intentional.
Slab buffers and MPI traffic use complex(c_float_complex)/MPI_COMPLEX while FFTs operate on double‑complex. This introduces implicit casts and potential precision loss. If double precision is required end‑to‑end, align buffer kinds and MPI datatypes.

🔧 Possible alignment (if double precision is desired)

-    complex(c_float_complex), allocatable :: data_cmplx_slabz(:, :, :), data_cmplx_slaby(:, :, :)
+    complex(c_double_complex), allocatable :: data_cmplx_slabz(:, :, :), data_cmplx_slaby(:, :, :)

-        call MPI_Alltoall(sendbuf_sf, NxC*Nyloc*Nzloc, MPI_COMPLEX, &
-                          recvbuf_sf, NxC*Nyloc*Nzloc, MPI_COMPLEX, MPI_COMM_WORLD, ierr)
+        call MPI_Alltoall(sendbuf_sf, NxC*Nyloc*Nzloc, MPI_DOUBLE_COMPLEX, &
+                          recvbuf_sf, NxC*Nyloc*Nzloc, MPI_DOUBLE_COMPLEX, MPI_COMM_WORLD, ierr)

Also applies to: 895-897, 933-935, 973-975, 1015-1017

src/pre_process/m_data_output.fpp (1)

842-865: Fix MPI datatype/count for levelset and levelset_norm writes.
Using mpi_p and raw element counts bypasses mixed-precision sizing; use mpi_io_p and scale by mpi_io_type.

🐛 Proposed fix

-                call MPI_FILE_SET_VIEW(ifile, disp, mpi_p, MPI_IO_levelset_DATA%view, &
+                call MPI_FILE_SET_VIEW(ifile, disp, mpi_io_p, MPI_IO_levelset_DATA%view, &
                                        'native', mpi_info_int, ierr)
-                call MPI_FILE_WRITE_ALL(ifile, MPI_IO_levelset_DATA%var%sf, data_size*num_ibs, &
-                                        mpi_p, status, ierr)
+                call MPI_FILE_WRITE_ALL(ifile, MPI_IO_levelset_DATA%var%sf, data_size*num_ibs*mpi_io_type, &
+                                        mpi_io_p, status, ierr)
…
-                call MPI_FILE_SET_VIEW(ifile, disp, mpi_p, MPI_IO_levelsetnorm_DATA%view, &
+                call MPI_FILE_SET_VIEW(ifile, disp, mpi_io_p, MPI_IO_levelsetnorm_DATA%view, &
                                        'native', mpi_info_int, ierr)
-                call MPI_FILE_WRITE_ALL(ifile, MPI_IO_levelsetnorm_DATA%var%sf, data_size*num_ibs*3, &
-                                        mpi_p, status, ierr)
+                call MPI_FILE_WRITE_ALL(ifile, MPI_IO_levelsetnorm_DATA%var%sf, data_size*num_ibs*3*mpi_io_type, &
+                                        mpi_io_p, status, ierr)

src/post_process/m_global_parameters.fpp (1)

1066-1070: Deallocate q_filtered_wrt MPI_IO buffers before nulling.
The extra stats slots are allocated but only nullified here, which leaks memory. Based on learnings, please use the device-aware deallocation macro.

🐛 Proposed fix

             if (q_filtered_wrt) then
                 do i = sys_size + 1, sys_size + volume_filter_dt%stat_size
-                    MPI_IO_DATA%var(i)%sf => null()
+                    if (associated(MPI_IO_DATA%var(i)%sf)) then
+                        @:DEALLOCATE(MPI_IO_DATA%var(i)%sf)
+                    end if
+                    MPI_IO_DATA%var(i)%sf => null()
                 end do
             end if

toolchain/mfc/case_validator.py (2)

1111-1124: Remove trailing whitespace in the new validators (lint failure).
Pylint is failing on trailing whitespace in these added blocks; trimming resolves the pipeline error.

Also applies to: 1162-1200

---

1194-1195: Fix the even-dimension parity check.
The current (dim + 1) % 2 test rejects valid even sizes; check the dimension directly.

🐛 Proposed fix

-            self.prohibit((m_glb + 1) % 2 != 0 or (n_glb + 1) % 2 != 0 or (p_glb + 1) % 2 != 0,
+            self.prohibit(m_glb % 2 != 0 or n_glb % 2 != 0 or p_glb % 2 != 0,
                          "volume_filter_momentum_eqn requires global dimensions divisible by 2 in every direction")

src/common/m_mpi_common.fpp (1)

386-391: Avoid overwriting base MPI views when initializing filtered views.
The loop starts at sys_size, which can reinitialize the last base variable’s view. Start at sys_size + 1 instead.

🐛 Proposed fix

-        do i = sys_size, alt_sys
+        do i = sys_size + 1, alt_sys
             call MPI_TYPE_CREATE_SUBARRAY(num_dims, sizes_glb, sizes_loc, start_idx, &
                                           MPI_ORDER_FORTRAN, mpi_p, MPI_IO_DATA%view(i), ierr)
             call MPI_TYPE_COMMIT(MPI_IO_DATA%view(i), ierr)
         end do

src/common/m_ib_patches.fpp (1)

619-665: Avoid shared-state coordinate conversion inside the GPU loop.
s_convert_cylindrical_to_cartesian_coord writes module-level cart_y/cart_z, which is unsafe in parallel loops. Use the pure function with thread-local temporaries instead.

🐛 Suggested fix (thread‑local conversion)

-                        if (grid_geometry == 3) then
-                            call s_convert_cylindrical_to_cartesian_coord(y_cc(j), z_cc(k))
-                        else
-                            cart_y = y_cc(j)
-                            cart_z = z_cc(k)
-                        end if
+                        if (grid_geometry == 3) then
+                            tmp_cyl = (/0._wp, y_cc(j), z_cc(k)/)
+                            tmp_cart = f_convert_cyl_to_cart(tmp_cyl)
+                            cart_y = tmp_cart(2)
+                            cart_z = tmp_cart(3)
+                        else
+                            cart_y = y_cc(j)
+                            cart_z = z_cc(k)
+                        end if

src/simulation/m_data_output.fpp (2)

85-123: Serial I/O path does not forward filtered data arguments.

When parallel_io is false (Line 117-118), the call to s_write_serial_data_files does not pass the new filtered arguments (filtered_fluid_indicator_function, stat_q_cons_filtered, etc.). If q_filtered_wrt is enabled with serial I/O, the filtered data will not be written.

Either extend s_write_serial_data_files to accept and handle these arguments, or add a runtime check/abort when q_filtered_wrt is true and parallel_io is false.

---

801-814: Argument count exceeds the ≤6 guideline.

s_write_parallel_data_files now accepts 10 arguments. As per coding guidelines, consider introducing a derived-type params struct to bundle the filtered/stat fields and reduce the signature complexity.

src/simulation/m_ibm.fpp (2)

625-634: Off-by-one risk: loop accesses s_cc(index + 1) but bound allows index == bound.

The bound calculation at lines 627, 630, 633 sets bound = m/n/p + buff_size, but the while-loop at line 648-656 accesses s_cc(index + 1). When index == bound, s_cc(index + 1) reads past the allocated array bounds.

🐛 Suggested fix

                 if (dim == 1) then
                     s_cc => x_cc
-                    bound = m + buff_size
+                    bound = m + buff_size - 1
                 elseif (dim == 2) then
                     s_cc => y_cc
-                    bound = n + buff_size
+                    bound = n + buff_size - 1
                 else
                     s_cc => z_cc
-                    bound = p + buff_size
+                    bound = p + buff_size - 1
                 end if

---

651-655: Replace error stop with s_mpi_abort for MPI safety.

Using error stop can leave other MPI ranks hanging. As per coding guidelines, use s_mpi_abort for error termination.

🐛 Suggested fix

                         if (index < -buff_size .or. index > bound) then
-                            print *, "proc_rank=", proc_rank, "temp_loc=", temp_loc, " index=", index, "ib=", patch_id, "dim", dim, "dir", dir, "i, j, k", i, j, k
-                            print *, "Increase buff_size further in m_helper_basic (currently set to a minimum of 10)"
-                            error stop "Increase buff_size"
+                            call s_mpi_abort("Image point search exceeded buffer bounds. Increase buff_size (minimum 10).")
                         end if

src/simulation/m_additional_forcing.fpp (1)

170-172: Division by zero risk when local density is zero.

Line 172 divides by rho, which is computed by summing fluid densities. At phase interfaces or in regions where all fluids are absent, rho could be zero or near-zero, causing division by zero or numerical instability.

🐛 Suggested fix

                     ! f_E
-                    q_periodic_force(3)%sf(i, j, k) = (P_inf_ref*gammas(1) - phase_eps)*forcing_dt &
-                                                      + q_cons_vf(contxe + mom_f_idx)%sf(i, j, k)*q_periodic_force(2)%sf(i, j, k)/rho
+                    if (rho > 1.e-12_wp) then
+                        q_periodic_force(3)%sf(i, j, k) = (P_inf_ref*gammas(1) - phase_eps)*forcing_dt &
+                                                          + q_cons_vf(contxe + mom_f_idx)%sf(i, j, k)*q_periodic_force(2)%sf(i, j, k)/rho
+                    else
+                        q_periodic_force(3)%sf(i, j, k) = (P_inf_ref*gammas(1) - phase_eps)*forcing_dt
+                    end if

src/simulation/m_volume_filtering.fpp (4)

214-216: Misleading error message: references 'p' but checks 'Ny'.

The check validates mod(Ny, num_procs) but the error message incorrectly states the requirement is for p.

✏️ Suggested fix

         if (mod(Ny, num_procs) /= 0) then
-            call s_mpi_abort('Volume filtering requires p to be divisible by the number of ranks')
+            call s_mpi_abort('Volume filtering requires Ny to be divisible by the number of ranks')
         end if

---

585-598: Critical: Out-of-bounds array access in finite-difference stencil.

The loop iterates i = 0, m, j = 0, n, k = 0, p, but the stencil accesses i + l where l ∈ [-4, 4]. At boundary cells (e.g., i = 0), this accesses fluid_indicator_function%sf(-4, j, k), causing out-of-bounds memory access unless idwbuff provides sufficient ghost cells.

Verify that idwbuff bounds provide at least 4 ghost cells in each direction, or restrict the loop ranges:

🐛 Option 1: Restrict loop ranges (safer)

         $:GPU_PARALLEL_LOOP(collapse=3, copyin='[fd_coeffs]')
-        do i = 0, m
-            do j = 0, n
-                do k = 0, p
+        do i = 4, m - 4
+            do j = 4, n - 4
+                do k = 4, p - 4
                     $:GPU_LOOP(parallelism='[seq]')
                     do l = -4, 4

🔍 Option 2: Verify buffer allocation

#!/bin/bash
# Check idwbuff allocation bounds to verify sufficient ghost cells
rg -n 'idwbuff' --type-add 'fpp:*.fpp' --type=fpp -C3 | head -60

---

785-787: Potential division by zero on filtered density.

When q_cons_filtered(1)%sf(i, j, k) (filtered density) approaches zero inside solid regions, this division produces NaN/Inf values.

🐛 Proposed fix with epsilon floor

+                            real(wp), parameter :: eps_rho = 1.0e-12_wp
                             reynolds_stress(l, q)%sf(i, j, k) = reynolds_stress(l, q)%sf(i, j, k) &
-                                                                - (q_cons_filtered(momxb - 1 + l)%sf(i, j, k)*q_cons_filtered(momxb - 1 + q)%sf(i, j, k)/q_cons_filtered(1)%sf(i, j, k))
+                                                                - (q_cons_filtered(momxb - 1 + l)%sf(i, j, k)*q_cons_filtered(momxb - 1 + q)%sf(i, j, k) &
+                                                                   / max(q_cons_filtered(1)%sf(i, j, k), eps_rho))

---

1552-1560: Multiple divisions by filtered_fluid_indicator_function risk NaN in solid regions.

Lines 1556, 1605, 1655, and 1707 all divide by filtered_fluid_indicator_function%sf(i, j, k), which is zero inside particles. This systematically produces NaN/Inf values throughout the filtered fields.

Apply a consistent epsilon guard across all four locations:

🐛 Proposed pattern for all divisions

Define a module-level constant and apply uniformly:

real(wp), parameter :: eps_indicator = 1.0e-12_wp

Then replace each division:

-q_cons_filtered(l)%sf(i, j, k) = data_real_3D_slabz(i + 1, j + 1, k + 1)*fft_norm/filtered_fluid_indicator_function%sf(i, j, k)
+q_cons_filtered(l)%sf(i, j, k) = data_real_3D_slabz(i + 1, j + 1, k + 1)*fft_norm &
+    / max(filtered_fluid_indicator_function%sf(i, j, k), eps_indicator)

Apply the same pattern at lines 1605, 1655, and 1707.

src/pre_process/m_global_parameters.fpp (1)

949-954: Use @:ALLOCATE for the new levelset MPI buffers.
These allocations are in a .fpp module and should use the device-aware macro. As per coding guidelines, please prefer @:ALLOCATE here.

♻️ Suggested change

-            allocate (MPI_IO_levelset_DATA%var%sf(0:m, 0:n, 0:p, 1:num_ibs))
-            allocate (MPI_IO_levelsetnorm_DATA%var%sf(0:m, 0:n, 0:p, 1:num_ibs, 1:3))
+            @:ALLOCATE(MPI_IO_levelset_DATA%var%sf(0:m, 0:n, 0:p, 1:num_ibs), &
+                       MPI_IO_levelsetnorm_DATA%var%sf(0:m, 0:n, 0:p, 1:num_ibs, 1:3))

src/post_process/m_global_parameters.fpp (1)

872-878: Use @:ALLOCATE for q_filtered_wrt MPI_IO_DATA buffers.
These new allocations are in a .fpp module; using the device-aware macro keeps GPU builds consistent. As per coding guidelines, please prefer @:ALLOCATE here.

♻️ Suggested change

-                allocate (MPI_IO_DATA%var(i)%sf(0:m, 0:n, 0:p))
+                @:ALLOCATE(MPI_IO_DATA%var(i)%sf(0:m, 0:n, 0:p))

src/post_process/m_data_input.f90 (1)

671-680: Consider bundling MPI I/O params into a derived type.
This signature has 8 parameters; grouping the MPI offsets into a small derived type would align with the ≤6-argument guideline and simplify call sites. As per coding guidelines, ...

src/simulation/m_ibm.fpp (2)

713-723: Remove or enable commented-out code.

Lines 714-723 contain commented-out alternative logic for 1D subsection checks. If this code is no longer needed, remove it to reduce clutter. If it's intended for future use, consider adding a TODO comment explaining its purpose.

---

805-811: Remove or enable commented-out code.

Lines 806-811 contain the same commented-out 1D subsection logic pattern. This duplicates the commented code in s_find_num_ghost_points. Either remove both or document the intent.

src/simulation/m_volume_filtering.fpp (3)

339-343: Consider using newunit= for automatic unit assignment.

Hardcoded unit number 100 may conflict with other file operations in the codebase. Using newunit= provides safer automatic assignment.

♻️ Suggested improvement

+        integer, save :: particle_force_unit
         ! file for particle forces
         if (compute_particle_drag) then
             if (proc_rank == 0) then
-                open (unit=100, file='particle_force.bin', status='replace', form='unformatted', access='stream', action='write')
+                open (newunit=particle_force_unit, file='particle_force.bin', status='replace', form='unformatted', access='stream', action='write')
             end if
         end if

Then update the write and close statements to use particle_force_unit instead of 100.

---

1037-1765: Subroutine s_filter_batch exceeds size guidelines (~730 lines).

Per coding guidelines, subroutines should be ≤ 500 lines. This subroutine contains highly repetitive FFT transform patterns for different field types (conservative vars, pressure, Reynolds stress, effective viscosity, interphase momentum exchange).

Consider extracting reusable helper routines:

1. s_fft_forward_xy: Handles X FFT → reorder → Y FFT → store to batch slab
2. s_fft_z_convolve: Handles Z FFT → convolution → inverse Z FFT
3. s_fft_backward_yx: Handles Y inverse FFT → reorder → X inverse FFT → store result

This would reduce s_filter_batch to ~100-150 lines of orchestration logic and improve maintainability.

---

3-2015: File exceeds 1000-line guideline (~2015 lines).

The coding guidelines recommend files ≤ 1000 lines. Consider splitting into focused submodules:

• m_volume_filtering: Public API, orchestration, core filtering workflow
• m_volume_filtering_fft: FFT plans, forward/backward transforms, slab transposes
• m_volume_filtering_terms: Stress tensor, Reynolds stress, effective viscosity computations